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Abstract 

Scattering of femtosecond laser pulses on resonant transmission and reflection gratings 
made of dispersive (Drude metals) and dielectric materials is studied by a time-domain 
numerical algorithm for Maxwell's theory of linear passive (dispersive and absorbing) me- 
dia. The algorithm is based on the Hamiltonian formalism in the framework of which 
Maxwell's equations for passive media are shown to be equivalent to the first-order equa- 
tion, d^f/dt = 7i^, where H is a linear differential operator (Hamiltonian) acting on a 
multi-dimensional vector ^ built of the electromagnetic inductions and auxiliary mat- 
ter fields describing the medium response. The initial value problem is then solved by 
means of a modified time leapfrog method in combination with the Fourier pseudospectral 
method applied on a non-uniform grid that is constructed by a change of variables and 
designed to enhance the sampling efficiency near medium interfaces. The algorithm is 
shown to be highly accurate at relatively low computational costs. An excellent agree- 
ment with previous theoretical and experimental studies of the gratings is demonstrated 
by numerical simulations using our algorithm. In addition, our algorithm allows one to 
see real time dynamics of long leaving resonant excitations of electromagnetic fields in 
the gratings in the entire frequency range of the initial wide band wave packet as well as 
formation of the reflected and transmitted wave fronts. 
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1 Introduction 



The purpose of the present study is twofold. First, we test a novel time-domain algorithm 
for the Maxwell's theory of linear, passive (dispersive and absorbing) media. The algorithm 
is based on (i) the Hamiltonian formalism for evolution differential equations [1], on (ii) the 
time leapfrog scheme [2], and on (in) the Fourier pseudospectral method [3] in combination 
with a change of variables that enhances the spatial grid resolution in designated domains 
(in the vicinity of medium interfaces) and, thereby, prevents the loss of accuracy due to the 
aliasing problem of the Fourier transform, while keeping the total spatial grid size fixed [4] , [5] . 
Boundary conditions at medium interfaces are not fixed in the algorithm, but rather medium 
parameters are allowed to have spatial discontinuities so that the correct boundary conditions 
are enforced dynamically [6], similarly to the wave packet method for quantum mechanical 
systems with discontinuous potentials. Although, each of these three main ingredients of our 
algorithm have been used individually in various computational problems in electromagnetism 
and quantum mechanics, to our knowledge they have never been put together in applications 
to numerical simulations of propagation of a wideband electromagnetic pulse in passive media 
and its scattering on targets made of dispersive and absorbing materials. By combining these 
methods, we have obtained an efficient, true time-domain algorithm that is highly accurate, 
which is a known virtue of pseudospectral methods of solving partial differential equations [7] . 

An essential advantage of time-domain numerical methods is that one can see all the imme- 
diate effects the medium and targets have on the propagating wideband wave packet. Yet, a 
single simulation of scattering of a wideband wave packet is sufficient to determine some basic 
electromagnetic properties of a target, e.g., transmission and reflection coefficients, in the entire 
frequency range of the initial wave packet. 

It has been reported [8] that a periodic thin-film metallic grating (either with holes or one- 
dimensional slits) can transmit more light at certain wavelengths than the projected area of 
the holes (or slits) in the grating would suggest, while at other wavelengths transmission is 
almost fully blocked. There is an ongoing discussion about the mechanism of such anomalous 
transmission. Among suggested mechanisms are formation of dynamical diffraction resonances 
in periodic metallic structures [9] , surface plasmons whose resonances are enhanced by the array 
of holes in the grating [10], and, specific to the slit gratings, open Fabry-Perot resonant cavities 
[11]. Since all the approaches produce essentially identical predictions for the transmission and 
reflection coefficients (the far-field) of the slit grating, it is, perhaps, necessary to have a closer 
look at the details of the electromagnetic field dynamics in the vicinity of the grating where 
deviations of theoretical predictions of a particular mechanism from the actual exact solution 
of Maxwell's equations might occur (see, e.g., a discussion in a recent work [12]). This is our 
second motivation of the present study. We apply our algorithm to scattering of a wideband 
electromagnetic pulse on various transmission and reflection slit gratings (dispersive metallic 
and purely dielectric ones). Our time-domain algorithm allows one to observe in details (in real 
time) formation and decay of long-living resonant excitations of electromagnetic fields as well as 
formation of resonant transmission and reflection wave fronts. Thanks to the use of the Fourier 
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pseudospectral method in combination with nonuniform grids, an extremely high accuracy of 
simulations can be achieved in the entire simulation volume and time span at relatively low 
computational costs. We think that these virtues of our time-domain algorithm would be useful 
for numerical studies of electromagnetic properties of other nanostructured materials [13]. 

The paper is organized as follows. Section 2 is devoted to the Hamiltonian formalism 
applied to the Maxwell's theory for general passive linear media. Maxwell's equations for elec- 
tromagnetic fields and medium responses are shown to be equivalent to a first order evolution 
differential equation similar to the Schrodinger equation in quantum mechanics. The wave 
function is a multidimensional column whose components are electromagnetic fields and aux- 
iliary fields describing the medium response. The Hamiltonian operator is a linear differential 
operator acting on a Hilbert space spanned by square integrable wave functions. In Section 
3, an example of the multi-resonant Lorentz model, which is widely used to describe passive 
media, is considered in the framework of the Hamiltonian formalism. In Section 4, we establish 
a relation between the norm of wave functions and electromagnetic energy. The discussion is 
limited to the Lorentz model and non- dispersive dielectric media. Section 5 is devoted to a 
general description of our algorithm. In particular, we show how the action of the Hamiltonian 
on wave functions is defined on a grid by means of the fast Fourier transform method. We prove 
that the Gauss law is enforced in the grid representation at no extra computational costs in our 
algorithm. Then we discuss how the sampling efficiency of the fast Fourier method can be en- 
hanced in designated spatial regions (typically at medium interfaces) by changing variables on 
the grid. The time evolution is done by a modified leapfrog method applied to the Schrodinger 
equation. We show that the conventional leapfrog method leads to an unstable algorithm for 
media with absorption and propose a general method to modify the leapfrog scheme to obtain 
a conditionally stable algorithm. Finally, we give an explicit realization of our algorithm in 
the case of the multi-resonant Lorentz model. Section 6 contains a detailed description of the 
actual computational scheme used in our simulations of extraordinary transmission and reflec- 
tion grating. We analyze the stability of the scheme using our general approach developed in 
Section 5. Section 7 is devoted to our numerical results. We also compare them with previous 
theoretical and experimental studies. Section 8 contains a brief conclusion. 

2 Maxwell's theory for passive media in the Hamiltonian 
formalism 

Let E and H be electric and magnetic fields, respectively, D and B the corresponding in- 
ductions, and P and M the medium polarization and magnetization vectors. Boldface letters 
denote three- vector fields in M 3 . Propagation of an electromagnetic wave packet in passive 
linear media in the absence of external radiating sources is described by the following set of 
equations [6] 



H 4> F (t) , 




3 



tf(t) = i> F (t) + i> R {t) = ^ F (t) + f dr X (t - r)^(r) , (2.2) 

(S).*--(S). ^- (£).«•- Ux T). w 

where c is the speed of light in vacuum, the overdot denotes the partial derivative d/dt with 
respect to time t, the spatial argument of the fields, denoted below by r, is suppressed. For 
generic anisotropic media, the medium response function x(t) is regarded as a linear operator 
(matrix) acting on the components of ip F . For isotropic media, it is a scalar. The response 
function is also position dependent for non- homogeneous media. Let ^(O) = ip F (0) be an 
initial wave packet with finite energy (finite L 2 (R 3 ) norm). We are interested in a finite norm 
solution of the initial value problem for Maxwell equations (2.1) subject to the constraint (the 
Gauss law) 

V • B(t) = V • D(t) = . (2.4) 

The response function x(t) is usually deduced from a microscopic model of the medium in 
question [6]. Therefore it is natural to assume that X (t) is a fundamental solution of some 
linear evolution differential equation so that 

C t ^ R (t) = u P i) F {t) , (2.5) 

where C t is a linear differential operator (a polynomial in d/dt) and uj v describes a coupling 
between applied electromagnetic fields and matter. In general, u g is a position dependent 
matrix. Causality of the medium response requires that the Fourier transform X (w) of the 
response function should have poles only in the lower part of the complex plane of uo [6] so 
that the Fourier transform of (2.2) reads ip R (u) = X (u)ip F \ou) . By taking the Fourier transform 
of Eq. (2.5), one finds that C(u)ip R (u) = u p ip F (u) and, hence, C t = u~ 1 \x(id/dt)]~ 1 . If 
the response function is known from measurements in the frequency domain, components of 
[x(w)] _1 can always be approximated by a polynomial with sufficient accuracy in a frequency 
range of interest. 

Now the Hamiltonian formalism [1] can be applied to (2.5) to transform it to an equivalent 
system of first-order differential equations 

i(t)='H F M m + VMFV(t) , (2.6) 

where £ is a column of auxiliary fields which are linear combinations of the response field and 
its time derivatives, save for the one of the highest order. There exists a linear operator 1Z such 
that ip R (t) = H£(t). Its explicit form depends on the details of going over to the Hamiltonian 
formalism. One can, for instance, identify the first component of £ with ip R . In this case, 1Z 
projects the column £ onto its first component. The response function can be expressed through 
the operators and Vmf and 1Z by making use of the fundamental solution of Eq. (2.6) 

X (t) = e^Ue^VMF , (2.7) 
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where 9{t) is the Heaviside function. Equation (2.7) can be regarded as a condition on possible 
choices of the operators H M and Vmf and 71. 

Applying 7Z to (2.6) we find 

iP{t) = nnU{t) y (2.8) 

where the relation TZVmf = has been used. It is not hard to be convinced that the latter 
relation holds when the first component of £ coincides with ip R . Any other choice of £ can 
be obtained by a canonical transformation [1] which is a linear nonsingular transformation of 
auxiliary fields, 

^S M £, det<S M ^0. (2.9) 

According to (2.6) and (2.7), 71 — > KS M X , 7i M — > SmHmS^ 1 , an d Vmf — > S m Vmf and, 
therefore, the condition TZVmf = remains true in the new basis of auxiliary fields. Denoting 
Vfm = —TZ7~L M and substituting (2.2) and (2.8) into (2.1) Maxwell's evolution equations (2.1) 
can be written in the Hamiltonian form 

V{t)=n^ F {t) + v FM m • (2.io) 

Finally, the electromagnetic and auxiliary fields are unified into one column (wave function) so 
that Eqs. (2.10) and (2.6) can be represented as a single first-order evolution equation 

V F (t) = H F V F (t) , (2.11) 

CO •*'=(£ ni)- « 2i2 » 

The index F indicates that electromagnetic fields are used as independent electromagnetic 
degrees of freedom. We shall refer to (2.12) as to a field representation. Accordingly, an 
induction representation is obtained by the similarity transformation 



(^^=S^ F , S=(l ^, H^S^S- 1 . (2.13) 



The corresponding blocks of Ti 1 have the form 

Ht = Ho, Vmi = Vmf, (2.14) 

Vim = V FM + nH F M -HoK=-HoK , (2.15) 
Km = n F M -V MF n. (2.16) 

In what follows we denote wave functions by ^(t) with Q being the representation index, F 
or I. This completes construction of the Hamiltonian representation of Maxwell's theory for 
linear passive media. 

Boundary conditions at medium interfaces and possible scatterers (targets) are not im- 
posed on electromagnetic fields, but rather they are enforced dynamically by allowing medium 
parameters to be discontinuous functions. The fundamental solution of (2.11) 

V F (t) = e tnF V F (0) , t > , (2.17) 
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is well defined for discontinuous "potentials" Vmf and Vfm> for example, by means of the 
Kato- Trotter product formula used in the path integral representation of (2.17) as shown in 
[14]. 



3 The Lorentz model 

The Hamiltonian formalism for the Lorentz model has been used in [15] to develop a finite 
differencing algorithm to study an electromagnetic pulse propagation in Lorentz media. Here 
we derive an explicit form of the Hamiltonian for the Lorentz model which is used in Section 5 
for the stability analysis of our algorithm. The Lorentz model of a passive medium is based on 
the assumption that the medium magnetization is zero, M = 0, while the medium polarization 
is described by a set of decoupled second-order differential equations [6] 

N 

P a + 2 7a P a + ^P a = ^ 2 a E, P = ^P a . (3.1) 

a=l 

where uj a are resonant frequencies, 7 a are damping coefficients, and u pa are plasma frequencies. 
As has been pointed out, no boundary conditions are imposed on electromagnetic fields at 
medium and/or target interfaces. Instead, the coupling constants uo, pa = uj pa (r) are allowed 
to have discontinuities at medium interfaces, or, from the physical point of view, they remain 
smooth but change rapidly, \ w \'Vup\/ujp y> 1, at the interface, where \ w is a typical wave 
length of the incident wave packet. The initial value problem is solved in the space of square 
integrable wave functions. Initial conditions for the response field are P(t = 0) = P(t = 0) = 0. 

Consider 2N real vector fields, £ j = 1,2, ...,2N, such that 

Pa = (WpaM»)£ 2 o-l » ( 3 - 2 ) 
£20-1 = ^2a , £ 2 a = - 2 7a£ 2 a ~ ^2a-l + ^P«E . (3.3) 

Thus, the original system of second order equations has been converted into the first order 
system. The operator 1Z is defined by (3.2). After simple algebraic transformations, we infer 

M , (3.4) 



Vmf = — V^ M , (3.5) 



Hm — diag (Hmi, 7~Cm2i ' ' ' > ^mn) > ^Ma — (^_ u ^ 



(3.6) 



where diag indicates that the corresponding matrix is block-diagonal with blocks listed in the 
order from the upper left to lower right corners. Note that the matrices Vfmo, and Tip Ma act 
on a six-dimensional column £ a composed of two vectors £ 2a -i an< ^ £ 2 <r Therefore they should 
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be understood as composed of 3 x 3 blocks. Each block is obtained by multiplying the unit 
matrix by the number indicated in place of the block in (3.4) and (3.6). 

Another convenient way to introduce the Hamiltonian formalism is to use N complex vector 
fields C a which satisfy the first order differential equation 

C« = - H»E , Pa = ^ (Ca + C«) , (3.7) 

where A a = — 7 a +i^ a and v a = - 7^. This representation is defined only if 7 a < uo a (i.e., the 
attenuation is not high). From the numerical point of view, solving a decoupled system of iV first 
order differential equation and taking complex conjugation (denoted here by an over bar) is less 
expensive than solving the original system of second-order differential equations for the medium 
polarization. In general, there is always a freedom of choosing a new basis for the auxiliary field 
space (2.9). If the evolution operator exp(t7i Q ) is computed in one basis, it can be computed in 
another basis by a suitable similarity transformation. This is an important observation because 
the auxiliary field basis can be chosen in a way that facilitates computation of the evolution 
operator, e.g., to speed up simulations. For instance, in the complex representation (3.7), the 
matter Hamiltonian Ti^ is diagonal. The corresponding transformation of auxiliary fields is 
given by 

(t')=£te (38) 

To transform the whole system into this representation, the Hamiltonian 7i F is replaced by 
S~ 1 7i F S and the wave function ^> F by S^ F where S is block-diagonal with the unit matrix in 
the upper left (field) corner and with Sm in the lower right (matter) corner. 



4 Energy and the norm of state vectors 

Accuracy and convergence of a numerical algorithm is defined relative to some norm. Let us 
discuss the choice of a norm in the space spanned by wave functions ^/ Q (t). The discussion 
is limited to the Lorentz model and the case of a nonhomogeneous, nondispersive medium 
(dielectric) which are used in our numerical simulations. 



I The Lorentz model 

Consider a multi-resonant Lorentz model with no attenuation 7 a = 0. The field and matter 
evolution equations can be obtained from the variational principle for the action 



= J dtL = Jdtjdr 



-(E 2 

2 v 



B 2 ) + i£ 



+ P E 



(4.1) 
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where the polarization of the medium is expressed via matter fields as P = ^2 ou pa 'd a . The 
electromagnetic degrees of freedom are described by the vector and scalar potentials, respec- 
tively, A and (p. The fields are defined by E = — V<£> — A and B = V x A. The units are 
chosen in this Section so that c = 1. The least action principle for the scalar potential (p leads 
to the Gauss law, V • D = 0, for the vector potential A to the Maxwell's equation, D = V x B, 
and for the matter fields "d a to the medium polarization evolution equation (3.1) of the Lorentz 
model with no attenuation, 7 a = 0. The second Maxwell's equation and the Gauss law for 
the magnetic field follows from the relation B = V x A by taking its time derivative and 
divergence, respectively. 

The energy of the system coincides with the canonical Hamiltonian which is obtained by 
the Legendre transformation [1] of the Lagrangian L for the velocities A and i9 a . The canonical 
momenta are iv a = <5L/<5$ a = $ a and II = 5L/5A = — E — P = — D. Doing the Legendre 
transformation, we find the canonical Hamiltonian (energy) of the system 



E(t) 



jdv f^7r a -^ a + n-Aj -L = ±Jdr 



E 2 + B 2 + ]>> 2 + u«) 



(4.2) 



where the Gauss law V • D = has been used. The energy conservation follows directly from 
the Noether theorem [1] applied to the time translation symmetry of the action (4.1), E(t) = 0. 
Equation (4.2) becomes the conventional expression for the electromagnetic energy in a passive 
medium [6] when 7r a and # a are replaced by the corresponding solutions of the equations of 
motion with initial conditions ir a (t = 0) = "d a (t = 0) = 0. 

An important observation is that the Noether integral of motion (4.2) coincides with the 
norm squared of the corresponding state vector 



E 



dr^ F ^ F 



(4.3) 



where pi — S 1 . This follows from the fact that, if we identify £ 2a = ir a and £ 2o -i = w «^ai 
the Hamiltonian equations of motion, i) a = 5E/5ir a and 7r a = —8E/8'd a , coincide with (3.3) 
when 7 a = 0. Note that canonically conjugated electromagnetic variables are A and — D. 
Therefore, the coupling between the electromagnetic and matter degrees of freedom in the 
Hamiltonian equations of motion is generated by the term E 2 = (D — P) 2 in (4.2). Thus, in 
the absence of attenuation, the norm of the state vector is proportional to the wave packet 
electromagnetic energy and, hence, is conserved. 

The norm conservation also follows from anti-Hermiticity of the Hamiltonian 7i F ^ = —7i F if 
7 a = 0, while (4.3) establishes a relation between the electromagnetic energy and the norm. In 
the induction representation, the norm in the measure space, defined by the operator \i in (4.3), 
is also conserved by construction. Consequently, the Hamiltonian is anti-Hermitian relative to 
the measure space scalar product, H^fi = —jiH 1 . 

The norm (energy) conservation can be used to control numerical convergence, especially 
when the aliasing problem in the fast Fourier transform is present, e.g., when parameters of 
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the medium are discontinuous functions in space. In a properly designed algorithm the loss of 
energy (norm) due to attenuation should be controlled by the symmetric part of the Hamiltonian 
operator 

E(t) = -J2^J dr &(*) = \ (* F (0,Vf < , (4.4) 

a 

where = V F = (H F ^ + H F )/2 < (a negative semidefinite operator) which is, in this case, 
a diagonal matrix with nonpositive elements. 

II Nondispersive media 

If the medium in question does not have dispersion and absorption, the formalism is simplified. 
Let e = e(x) be the dielectric constant of the medium, D = eE. If the medium is not isotropic, 
then e is symmetric positive definite 3x3 matrix everywhere in space. We rewrite Maxwell's 
equations in the form 

fi(t) = Ua^it) , H G = ( _ cV ^ } ) , (4.5) 

where the parenthesis in (e~ l ) mean that the induction is first multiplied by e^ 1 and then the 
curl of the resulting vector field is computed. Consider the scalar product 

Wi^M =Jdr ^W 2 7 , K=( e Q J) ■ ( 4 - 6 ) 

The Hamiltonian is ant i- Hermit ian with respect to this measure space scalar product, 7i G fi £ = 
—fi £ Ti,G- Therefore the corresponding norm is preserved in the time evolution generated by 
exp(tHc), that is, (vp 1 (t) , (x^ 1 (t)) = (^(O), /i e ^ / (0)). The electromagnetic energy of the wave 
packet is conserved because it is proportional to the measure space norm of the initial state 
vector. 



5 The algorithm 

I The grid representation 

Consider an equidistantly spaced finite grid with periodic boundary conditions. Let Ar be the 
grid step and n be a vector with integer valued components. Then the dynamical variables 
are projected onto the grid by taking their values at grid points r = nAr, that is, the wave 
function is replaced by a vector (column) whose dimension is determined by the grid size, 
vj/Q(r) — > $^(nAr) = \f^. A cubic grid is assumed. Consider a discrete Fourier transformation 
associated with the grid, ^ Q (nA; ) = £n' -7w*2 = where = = 1 - The 
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reciprocal lattice step is k = 2ir/Ar. The grid spatial size L and step Ar must be chosen so that 
the Fourier transform of the initial wave packet has support within the region k e [k min , k max ] 
where k — |k|, k max = k and k min = 2n/L. The action is of any position dependent operator 
V Q = V Q( r j is defined by 

VQ(r)tf?(r) -> V Q (nAr)^°(nAr) = E V nn>^n> > 

n' 

where V^ n , = 5 nn /V Q (nAr) is a diagonal matrix. Let Hq = 7Y^(V) depend only on the V 
operator. The projection of its action onto the grid is then defined via the discrete Fourier 
transform 

W?(V)*?(r) - E^nn^S , «!L = E (^W • (5.2) 

r=nAr ' ^- ~ ' 

n' n" 

The projection (5.2) is performed by the fast Fourier transform method. 

In what follows, the action of a product of and Ti^ on any state vector is understood 
as multiplication of ^ by the corresponding matrices, defined in (5.1) and (5.2), in the order 
specified in the product. The main advantage of using the Fourier basis is the exponential con- 
vergence (versus the polynomial one in finite differencing schemes) [7] as the grid size increases, 
which allows one to substantially increase the accuracy of simulations. 



II The Gauss law 

Another advantage of the Fourier basis is that the Gauss law is enforced at no extra compu- 
tational cost. In the grid representation defined above, the Gauss law (2.4) requires that the 
Fourier transforms of the inductions D(k) and B(k) remains perpendicular to the reciprocal 
grid vector k = n/c at any moment of time. In our algorithm, as we shall show shortly, the 
time evolution is generated by applying powers of the Hamiltonian to the wave function. In 
the induction representation, the action of powers of H 1 always produces the cross product 
k x C(k), for some C(k) regular at k = 0, in the entries of ^ T that correspond to the elec- 
tromagnetic inductions. Hence, in the grid representation the wave function (H 1 )" 1 ^ 1 satisfies 
the Gauss law for any power m because of the trivial identity k • (k x C(k)) = valid for any 
vector k of the reciprocal grid. 

It should be noted that a high accuracy of the Gauss law is essential to achieve a high 
accuracy of simulated electromagnetic fields near medium interfaces. 



Ill Improving sampling efficiency by changing variables 

As is well known from the Fourier analysis, the convergence rate can be affected for functions 
which have discontinuities [3]. The latter is, unfortunately, the case in electromagnetic scat- 
tering problems [6]. Suppose there is an interface between two media. It can be deduced from 
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the dynamical Maxwell's equations that the components of electric and magnetic fields, E and 
H, tangential to the interface must be continuous, provided there is no surface electric current 
on the interface. From the Gauss law it follows that components of the inductions, D and B, 
normal to the interface must be continuous, provided there is no surface charge on the inter- 
face. In contrast, normal components of the fields and tangential components of the inductions 
can be discontinuous. Their discontinuities are determined by discontinuities of medium pa- 
rameters (e.g., discontinuities in plasma frequencies in Lorentz models). Therefore, in either 
the induction or field representation, there are components which suffer discontinuities at the 
interface. The only way to cope with the problem, while keeping the use of the Fourier basis, 
is to make the grid finer [3]. This would lead to a substantial waste of computational resources 
because the conventional fast Fourier transform is defined on a uniform periodic grid, while the 
sampling efficiency should only be enhanced in the neighborhood of medium interfaces. The 
use of wavelet bases might be helpful for such a task [16] in time domain algorithms. Here we 
retain the Fourier basis, and increase the sampling efficiency by a change of variables. 

The basic idea can be understood with a one- dimensional example. Let z be a physical co- 
ordinate. Consider a change of variables defined by z = f(y) where y is an auxiliary coordinate. 
An equidistant grid y n = nAy, with n being integers, of the auxiliary coordinate generates a 
non-uniform grid of the the physical coordinate, z n = f(nAy). Assuming Ay to be sufficiently 
small and f(y) sufficiently smooth, the physical grid spacing can be approximated as 

Az n = z n+1 - z n w Ay f (nAy) . 

So, if f'(y) = 1, then Az n = Ay and the grid is equidistant. By making the derivative 
< f'(y) < 1 in some designated areas, one can achieve a desired local grid density in the 
physical space, while keeping the total grid size fixed. For example, if it is necessary to increase 
the sampling efficiency in the vicinity z = 0, one can take f'(y) = 1 — a [l + b^y 2 ]^ 1 , where 
< a < 1, and , hence, 

f(y) = V - ^tan -1 (6y) = y-g(y,a,b) . 

By adjusting parameters a and b, the local grid density can be changed as desired. Consequently, 
if the sampling efficiency is to be enhanced at several points yi, a suitable change of variables 
can be of the form f(y) = y — g(y — y i: bi). The fast Fourier transform is applied on a 
uniform grid of the auxiliary variable. 

Upon the change of variables, the integration measure in the scalar product and the deriva- 
tive transform, respectively, dz = dyf\y) and d z = [/'(y)] -1 ^. When projected on a uniform 
grid of the physical coordinate by means of (5.2), the derivative operator d z becomes an anti- 
Hermitian matrix. When the rule (5.2) is applied on a uniform grid of the auxiliary coordinate 
y, the derivative operator d z is no longer represented by an anti-Hermitian matrix, although it 
is still an anti-Hermitian linear operator, but in the measure space where the scalar product 
is defined with the weight f'inAy) at each lattice cite. From the numerical point of view it 
is convenient to have an explicitly anti-Hermitian matrix representation of the derivative on 
the grid. Due to round-off errors, the exact anti-Hermiticity of d z in the measure space can 
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be violated in the grid representation, which, in turn, may lead to numerical instabilities of 
simulations. Note that if a computed matrix is known to be anti-Hermitian, then only half of 
its elements is to be computed, while the other elements are restored by the symmetry. In the 
measure space, an explicit anti-Hermiticity of a linear operator is much more difficult to main- 
tain in numerical simulations because the symmetry relation between matrix elements depends 
on the scalar product measure. For this reason, the wave function is rescaled by the square 
root of the Jacobian [5], \P — > vT 7 ^ so that the integration measure becomes dy leading to the 
conventional Euclidean scalar product in the grid representation (with the uniform unit weight 
at each grid site). In such a representation, the derivative operator d z — > [/'J -1 / 2 ^/'] -1 / 2 
becomes again an explicitly anti-Hermitian matrix in the grid representation defined by the 
rules (5.1) and (5.2). 

IV A modified temporal leapfrog scheme 

A conventional temporal leapfrog method of solving the Schrodinger equation is based on the 
iteration scheme [2] 

^ Q (t + At) = ^ Q (t - At) + 2AtH Q ^ Q (t) , (5.3) 

so that the wave function in the consecutive time moment is computed from the wave function 
at two previous moments of time, where At is the time step. The action of the Hamiltonian 
is computed by pseudospectral methods, in particular, by means of the Fourier basis and 
the fast Fourier transform in our algorithm. The leapfrog algorithm is conditionally stable for 
an anti-Hermitian H Q , which is true for nondispersive media, but not the case for media with 
absorption. 

To investigate the stability, let us introduce the amplification matrix for the leapfrog algo- 
rithm, ^(t + At) = G(At)^(t) (the representation index, Q, is suppressed for a moment). From 
(5.3) it follows that Q{M) satisfies the equation G 2 (At) - 2AtHG{At) -1=0 which has two 
solutions 

(At) = HAt ± Vl + H 2 At 2 . (5.4) 

The stability of the algorithm requires that the energy norm of both approximate solutions 
[Q^ (At)] n ^ (0) must be uniformly bounded for n > 0. The necessary condition (but not 
sufficient) is the von Neumann condition that the spectral radius of the amplification matrix 
does not exceed 1. If a complex number Re t(fi is an eigenvalue of AtTi, then the von Neumann 
condition implies that tp = ±7r/2 and R 2 < 1. In other words, eigenvalues of AtTi must 
be imaginary and their magnitude should not exceed 1. If in addition we demand that the 
Hamiltonian is diagonalizable, then a conditional stability can be achieved for sufficiently small 
At. Indeed, in this case there exists an non-singular S such that 

H = S^HsS , H\ = -U s , (5.5) 

and AtHs satisfies the von Neumann condition. Since H and TCs have the same eigenvalues, the 
amplification matrices for H and TCs are related by the same similarity transformation (5.5). 
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Hence, the norm of a wave function obtained by the action of powers of (5.4) on an initial wave 
function is uniformly bounded. Note also that the Hamiltonian (5.5) is ant i- Hermit ian relative 
to the measure space scalar product, H)\i = —fJ>H, where \i = S'^S' 1 . The stability can also 
be proved via the equivalence of the conventional Euclidean norm and the /i-norm [14]. 

In the case of nondispersive media, the Hamiltonian Hg is ant i- Hermit ian in the measure 
space, and the von Neumann condition is fulfilled if 

Atck s max <l, (5.6) 

where h £ max the maximal norm of all wave vectors in the medium which can be estimated by 
y ' p{e)k max with k max being the maximal wave vector of the initial pulse in vacuum and p{e) 
the maximal spectral radius of the symmetric matrix e(x) over x (or simply the maximum of 
e(x) if the medium is isotropic). This can be understood from the following principle [17]. A 
finite difference scheme with variable coefficients is stable if all the corresponding schemes with 
frozen (i.e., fixed to a particular value everywhere in space) coefficients are stable. 

The von Neumann condition cannot be met if absorption is present because eigenvalues 
of the Hamiltonian must have real parts in order to account for exponential attenuation of 
field amplitudes. To circumvent this difficulty, the leapfrog scheme is modified in the following 
way [14]. We assume the Hamiltonian to be diagonalizable. The lack of eigenvectors of the 
Hamiltonian typically leads to solutions whose amplitudes grow polynomially in time [14]. This 
feature cannot be present in a physically reasonable model of passive media. So our assumption 
is justified from the physical point of view and, yet, the Lorentz model Hamiltonian is indeed 
diagonalizable. Let 

>H = H + V, (^,V^)<0 (5.7) 

for any \I/ in the Hilbert space and the von Neumann condition is satisfied for Ho, i.e., Ho has 
imaginary eigenvalues. The split (5.7) can be achieved in many ways. For instance, Ho can be 
obtained from H by setting all parameters responsible for attenuation to zero. In a physically 
acceptable model, the energy must be conservative in an absorption free medium and so must 
be the energy norm of wave vectors, and, therefore, the corresponding Hamiltonian must be 
ant i- Hermit ian (relative to the energy induced scalar product). In the Lorentz model this is 
easily seen in the field representation if we set Ho = H\ la= o which is explicitly anti-Hermitian 
(see Section 3). Then V is diagonal with matrix elements being zeros and — 2^y a . Another 
possibility is to identify Ho with the Hamiltonian in the vacuum, then V = H — Ho must be 
negative semidefinite in order to model exponential attenuation in passive media. Finally, one 
can also split the Hamiltonian into the sum of Hermitian and anti-Hermitian parts. 

After choosing a suitable split (5.7) we make a substitution \P(t) = exp(tV)<I ) (t) in the 
original evolution equation (2.11). The new wave function $(t) satisfies the equation with a 
time dependent Hamiltonian 

$(t) = e - tv Hoe tv <$>(t) = H(t)$(t) , (5.8) 

to be solved with the same initial condition $(0) = ^(0). Applying the leapfrog method to 
(5.8) we get $(t + At) = $(i - At) + 2AtH(t)$(t) valid up to 0(At 3 ). Returning to the initial 
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variables, we arrive at the following recurrence relation 

#(t + At) = £(2At)#(t - At) + 2AtC(At)H ^(t) , (5.9) 

where C(At) = exp(AtV). The amplification matrix, V(t + At) = Q c (At)^(t), for the recur- 
rence (5.9) satisfies the equation 

Gc(At) = C(2At)g c 1 (At) + 2At£(At)H ■ (5.10) 

A deviation of the approximate solution ^2(At) x I / (0) from the exact solution relative to the 
energy norm is of order At 2 for any n > 0. Thus, the scheme is convergent and, hence, a 
conditional stability exists for a sufficiently small At > according to a general theorem of 
Kantorovich [2] that establishes a general equivalence between convergence and conditional 
stability. The conditional stability of (5.9) can be understood from the following observation. 
Solving (5.10) by perturbation theory in At, it is not hard to find that 

G c (At) - Q v {At) = At 3 /C(At) , G v (At) = £(At/2)£ (At)£(At/2) , (5.11) 

where JC(At) is regular in the vicinity of At = and vanishes whenever 7i and V commute, 
and Q (Ai) is the amplification matrix when V is set to zero. The von Neumann condition is 
satisfied for Q Q {At) for a sufficiently small At > 0. Hence, powers of Q Q {At) applied to ^(0) 
cannot produce any exponential norm growth. Powers of Qv(At) differ from those of Qo(At) 
by factors that are powers of e A<v and, hence, can only produce exponential attenuation of the 
amplitude. Indeed, let V v (t) = e* v ^(0). Then d/dt(^ v ,^ v ) = 2(^ V ,V^ V ) < since V is 
negative semidefinite. Thus, the approximate solution produced by the amplification matrix 
Gv(At) has no exponential growth, while differing, relative to the energy norm, from that 
produced by Qc{At) by order of 0(At 2 ). Therefore the modified leapfrog scheme can be made 
conditionally stable and as accurate as desired by reducing the time step. 

It should be noted that our arguments do not prove that there cannot be any exponential 
growth of the norm ||\l/(£)|| = (^(t), ^(t)) 1 / 2 in the modified leapfrog scheme (5.9). All we 
can claim is that ||^(t)|| < exp(Kt) which is also true for the conventional leapfrog scheme. 
The difference is that in the modified leapfrog scheme K ~ 0(At 2 ) as we have argued (a 
consequence of (5.11) and uniform boundedness of powers of Qy\ while in the conventional 
leapfrog scheme the constant K is independent of At. Hence a possible exponential growth 
cannot be suppressed by reducing the time step in (5.3), while it can be done in (5.9). 



V An example of the Lorentz model 

To illustrate our general method we give an example of the Lorentz model commonly used to 
describe dispersive media. In the field representation of the Hamiltonian for the Lorentz model, 
we make the following decomposition 

HF =(Z VF o") + (o 4) sWf + VF - (512) 
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Substituting this decomposition into (5.9) we arrive at the following scheme 



ij F (t + At) = ij F (t - At) + 2AtH ip F (t) + 2At VFMaCit) , (5.13) 

a 

Cit + At) = e 2Am * C(t - At) + 2At e Atn M<> V MFa V(t) , (5.14) 



g'^Ma — g T»* 



_ sinh z/ a t , F x 

COsh I/ * + z y^-Ma + la) 



(5.15) 



where z> a = (7^ — u 2 ) 1 / 2 and the six-dimensional columns £ a are introduced in Section 3. The 
exponential (5.15) is easy to compute by expanding 'H F ia in the Pauli matrix basis (a basis for 
the Lie algebra su(2)) and then by using the well known formula for the exponential of a linear 
combination of Pauli matrices. For small attenuation, 7 a < cu a , we get u a = iu a . The hyperbolic 
functions in (5.15) become the trigonometric ones and u a is replaced by u a . Eigenvalues of the 
matter Hamiltonian are A a = — 7 a ± v a . Hence, ReA a < and amplitudes of the matter fields 
are always exponentially attenuated as t — > 00, unless 7 a = leading to ReA a = 0. 

The stability is ensured if Ti F satisfies the von Neumann condition. Let k max be the maximal 
norm of all wave vectors of the initial wave packet and u)™ ax be the maximal value of u p = 
(Z~2 a ^pa) 1 ^ 2 as a function of position, then a sufficient condition for stability reads 



^t\jc 2 kl ax + < 1 . (5.16) 

Here the idea of the frozen coefficients [17] has been used again. The left hand side of inequality 
(5.16) is nothing but the spectral radius of At7i F with frozen plasma frequencies so that 
uj p = oj™^. Note that it is not difficult to solve the characteristic equation for 7i F with frozen 
plasma frequencies by using the Fourier basis. The scheme (5.9) becomes especially simple in 
the case of small attenuation 7 a < uo a . In the complex representation of the auxiliary fields 
(3.8) (cf. (3.7)) the matter Hamiltonians H F Ia are diagonal and the action of its exponential is 
reduced to multiplication by a complex number e lUaAt . 

Finally, it should be mentioned that, by rearranging operators in the split, namely, by 
moving Vfm and Vmf to V F in (5.12), the stability condition (5.16) can be weakened to 
Atck max < 1. This would come at the price of having a more complicated expression for C(At). 
In the case of the Lorentz model it can still be computed analytically. The new split can also be 
viewed as the use of the induction representation in the modified leapfrog scheme, Ti 1 = Ti^ + V 1 
where Ti^ contains only the blocks of Ti 1 with the V operator. The proof of the weaker stability 
condition can be found in [14]. 



6 Metal gratings in the Drude formalism 

Here we apply our method to gratings made of a metal whose optical properties are described 
by the Drude formalism. This is an actual numerical scheme used in simulations in Section 7. 
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The metal dielectric constant as a function of frequency is given by 



eM = 1 + xH 



cot, 



uj{uj + irj) 



where u p is the plasma frequency, which is zero in the vacuum and constant in the region 
occupied by the metal as shown in Fig.l, and 77 is the absorption. The model coincides with 
a one-resonant Lorentz model if u = and 77 = 27. To satisfy the Gauss law exactly in 
simulations, we use the induction representation according to Section 5. II. An auxiliary field 
is chosen so that its first component equals P and the second is denoted The Hamiltonian 
evolution equations are taken in the form 



D = cV x B , 

B = -cV x (D - P) , 

P = Tl£, 



Z = -fit--*- (D-P) . 
To apply the modified leapfrog scheme the Hamiltonian is split into the sum 



/ 





cVx 





\ 




-cVx 





cVx 
















-77 


V 


2 —i 
-up] 







-77/ 



Ti 1 



V 1 = diag(0, 0, 0, -77) . 



K + V 1 , 



(6-1) 



Clearly, V 1 is negative semidefinite because 77 > 0. The stability of the modified leapfrog scheme 
requires that eigenvalues of Hq have zero real parts. This is indeed the case. In the Fourier 
basis, the V operator becomes ih. It is straightforward to find the characteristic polynomial 
det(7io — A). Its nonzero roots are A = ±7iy/c 2 k 2 + ~iJl. Hence the scheme is stable if the time 
step is chosen so that the condition (5.16) is satisfied. In particular, k max can be set to k being 
the step of the reciprocal lattice and u™ ax is the plasma frequency of the metal (silver in our 
simulations, see Section 7). An explicit form of the modified leapfrog scheme for the split (6.1) 
reads 



D(t + At) = D(f - At) + 2cAtV x B(t) , 

B(t + At) = B(t - At) - 2cAtV x [D(t) - P(t)] , 

P(t + At) = P(t - At) - 2r)At£(t) , 

£(t + At) = e' 2r,M t(t- At) - 2Atiu 2 p r ] - 1 e'^ [D(t) - P(t)] . 

The action of the curl is computed by the fast Fourier method in combination with a change of 
variables that enhances the sampling efficiency near the metal-vacuum interface. The details 
are in Section 7. 
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7 Results for extraordinary transmission gratings 



To test our method we applied it to study transmission properties of metal and dielectric 
gratings suspended in vacuum. Transmission and reflection gratings have been subject of a 
number of experimental and theoretical works [9]— [12], [18]. The interest is stimulated by nearly 
100% transmission or reflection within narrow wavelength range with possibility of using such 
grating as efficient filters. Moreover, extraordinary optical transmission has been observed in 
the 2D hole arrays [8], further stimulating theoretical and experimental interest to transmission 
properties of nanostructured materials [13]. 

In Fig.l we schematically represent our system comprising a metallic or dielectric slab with 
gratings. D is a grating period, a the size of the grating and h its thickness. For the sake of 
comparison with previously published results we have chosen D = 1.75/im and a = 0.30/xm. 
Transmission and reflection coefficients are computed as a function of the grating thickness 
h. An incident electromagnetic wave packet propagates along the ^-direction, normal to the 
slab. The polarization is such that the electric field vector is perpendicular to the gratings, 
while the magnetic field is parallel to them (the so called p -polarization). Calculations are 
performed in a finite (x,z) box of the size -15D < z < 12D, —D/2 < x < D/2. A uniform 
mesh of typically 256 knots is used in the x-direction and a nonuniform mesh of 512 knots 
generated by the change of variables is used in the ^-direction. The change of variables is 
used to enhance the sampling efficiency near the two interfaces in the z-direction. Periodic 
boundary conditions are insured in x through the pseudospectral approach based on the Fast 
Fourier Transform. To suppress an artificial reflection of the wave packet, absorbing layers 
are introduced at the box boundaries z = ±15D. The initial wave packet is Gaussian, which 
allows us to obtain transmission (reflection) coefficients within the entire frequency bandwidth 
of the initial wave packet by a single simulation. In what follows we are mainly interested in 
transmission (reflection) of the radiation with wavelengths larger than the grating period (zero 
order diffraction). Thus, reflected or transmitted waves propagate in the direction of the z-axis, 
the same as the incident radiation. Note, however, that this is not a limitation of our method 
which allows for a priori extraction of the entire scattering matrix for all wave vectors. 

In Fig. 2 we show transmission and reflection coefficients obtained for metallic gratings of 
variable depth. The dielectric response of the metal is described within the Drude formalism 
(see Section 6). Following [10] we use u p = 9eV and rj = O.leV representative for silver. As 
clearly seen in the Figure, the transmission coefficients exhibit narrow resonances for certain 
wavelengths. With increasing thickness of the gratings, h, the number of resonant structures 
increases. Our results are in a full agreement with previously published theoretical studies [10]. 
The only difference being, that the model used in [10] assumes perfect metal surfaces inside 
the gratings, neglecting possible absorption of the radiation. This leads to 100% transmission 
at the resonant frequencies. In our case, a part of radiation is absorbed by the metal so that 
the transmission never reaches 100%. 

The observed structures in transmission are associated with resonant modes of the elec- 
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tromagnetic field produced by coupled surface electromagnetic modes (called surface plasmon 
polaritons) and waveguide modes inside the gratings. Some of these resonances posses rela- 
tively long lifetimes. This can be immediately inferred from their width, as, e.g., the resonance 
located at A = 1.1D in the case of the grating with thickness h = lA/j,m. Another way to 
observe the trapped field resonances is to look at the time dependence of the field transmitted 
in the z-direction. The signal in Fig. 3 is registered by a detector placed at the distance of 
3.5-D behind the gratings. As seen in the figure, as soon as the resonances are populated by 
an incident 25/s-long pulse, they radiate the field during at least 125/s. Here the radiation 
time is determined by the lifetime of the narrowest resonance located at A = 1.1D. It is worth 
mentioning that, as shows the sum of the transmission and reflection coefficients, the total 
absorption is largest at the resonance positions, i.e., when the interaction time between the 
radiation (trapped mode) and the metal is large. 

Finally, Fig. 4 shows the field structure of a trapped mode corresponding to the narrowest 
resonance observed with the h = 2A/j,m thickness grating. The E x component of the electric 
field is presented. It is obtained via sufficiently long time propagation so that contributions 
from less long-lived states vanish. 

In Figs. 5 and 6 we show reflection coefficients for dielectric gratings suspended in vacuum. 
Simulations are performed using the conventional leapfrog scheme (5.3) applied to (4.5). Since 
the attenuation is absent, the scheme is stable if the time step satisfies (5.6). The dielectric 
material is modeled through the frequency independent dielectric constant e = 2 (Fig. 5) and 
e = 4 (Fig.6). Only the reflection coefficient is shown here since there is no absorption of the 
radiation in dielectric so that the transmission can simply be inferred from the unitarity of the 
scattering matrix. Without gratings the dielectric slabs are basically transparent in both cases. 
Introducing grating structures results in a complete reflection of the incident radiation within 
an extremely narrow wavelength bandwidth. The associated guided mode resonances have been 
extensively discussed in the literature [18]. In the case of e = 2 the resonances are so narrow 
that extraction of frequency dependent transmission and reflection coefficients by the Fourier 
transform of the scattered wave would require too large propagation time (see also Fig. 7). We 
had to stop our wave packet propagation before the radiation emitted by resonances ended. 
This explains why the reflection coefficients do not reach its maximal value 1 in this case. 
Consistently with the metal case, the number of resonant structures increases with increase 
of the width h of the gratings. The width of resonances increases with e as follows from the 
comparison of Figs. 5 and 6. Note also that resonances are associated with Fano profiles that 
usually arise because of the interference between the non-resonant and resonant contributions 
to the scattered wave. Such narrow reflection structures in the case of dielectric gratings have 
been usually studied by stationary methods in the frequency domain. An important advantage 
of the time-dependent study is that one has an immediate access to all the details the temporal 
evolution of electromagnetic fields in any desired part of the system. 

In Fig. 7 we show the time dependence of the field transmitted in the z-direction for gratings 
characterized by e = 2 and h = 0.8/im. With these parameters there is only one resonance in 
the reflection spectrum. The signal is registered by a detector placed at the distance of 3.5.D 
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behind the gratings. First, we observe that a 25/s incident pulse is transmitted through the 
structure without modification. The lasing effect, when the transmitted field is followed by 
basically monochromatic radiation, can clearly be seen. For readability reasons we could not 
show the complete time evolution in the figure, but the lasing effect lasts for at least (!) 2ps 
reflecting an extraordinary long lifetime of the trapped resonant field. It is this radiation which 
comes in the phase opposite to the corresponding harmonic in the transmitted initial pulse and 
leads, finally, to the zero transmission at the corresponding frequency. The same lasing effect to 
the left from the grating structure is responsible for the 100% reflection at the same frequency. 

In Fig. 8 we show a typical structure of the field corresponding to a trapped (resonant) 
mode. The D x component of the electric induction is represented in the case of dielectric 
gratings with e — 4. The thickness of the dielectric slab is h — 0.6/im. Note that in contrast to 
the metal grating structure, the field in the present case occupies the entire slab and not only 
the vacuum part of the grating. 



8 Conclusions 



We have developed a time domain algorithm for the initial value problem for the Maxwell's 
theory of linear passive media. The algorithm is based on (i) the Hamiltonian formalism for 
evolution differential equations, (ii) the Fourier pseudospectral method in which the sampling 
efficiency in designated space regions is enhanced by a suitable change of variables, and (Hi) 
the modified leapfrog scheme. We have analyzed the stability of the algorithm and found 
explicit stability conditions when passive media are described by multi-resonant Lorentz models. 
We have implemented and tested our algorithm for extraordinary transmission and reflection 
gratings whose optical properties have been studied in a number of theoretical and experimental 
works. Numerical simulations based on our algorithm are shown to produce extremely accurate 
data for the well studied far-field (zero-order diffraction) at relatively low computational costs. 
A single simulation of an incident wideband wave packet is sufficient to determine transmission 
and reflection properties of the gratings in the frequency range of the initial wave packet. In 
addition, our algorithm allows us to see a real time dynamics of formation of long-living resonant 
excitations of electromagnetic fields in the grating as well as formation of transmitted and/or 
reflected wave fronts in the entire frequency range of the initial wave packet. It is believed that 
our algorithm would be useful for numerical studies of other nanostructured materials. 
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Figures 



Fig. 1 A schematic representation of the studied system. The incident wave packet propa- 
gates along the normal to the slab containing the gratings (z-direction). 

Fig. 2. Calculated zero-order transmission (solid lines) and reflection (dashed lines) coef- 
ficients for metallic gratings described in the text. The results are presented as a function of 
the wavelength A of the incident radiation measured in the units of the period of the gratings, 
D. Different panels of the figure correspond to the different thickness h of the gratings, as 
indicated. 

Fig. 3. The electric field measured by a detector placed behind the metallic gratings 
with thickness h = lAfiin. Only the field corresponding to the zero-order transmitted wave 
propagating along the z-axis is represented. It is obtained by the Fourier analysis of the x- 
coordinate dependence of the field at the detector position. The signal is shown as a function 
of time measured in femtoseconds. 

Fig. 4. A snapshot of the x-component of the electric field, E x for the grating thickness 
h = 2.4/im. The results are presented as a function of x and z coordinates measured in units 
of D. Red and blue colors correspond, respectively, to positive and negative values of the field. 
The snapshot has been produced after a sufficiently long propagation time so that the field 
pictured in it indeed corresponds to the long-lived resonance giving enhanced transmission at 
A = 1.15D (see Fig.2). 

Fig. 5. Calculated zero-order reflection coefficient for dielectric gratings with e — 2. The 
results are presented as a function of the wavelength A of the incident radiation measured in the 
units of the period of the gratings, D. Different panels of the figure correspond to the different 
thickness h of the gratings, as indicated. 

Fig. 6. Calculated zero-order reflection coefficient for dielectric gratings with e — 4. Results 
are presented as a function of the wavelength A of the incident radiation measured in the units 
of the period of the gratings, D. Different panels of the figure correspond to the different 
thickness h of the gratings, as indicated. 

Fig. 7. The electric field measured by a detector placed behind the dielectric grating 
with thickness h = 0.8/im and dielectric constant e — 4. Only the field corresponding to the 
zero-order transmitted wave propagating along the z-axis is represented. It is obtained by the 
Fourier analysis of the ^-coordinate dependence of the field at the detector position. The signal 
is shown as a function of time measured in femtoseconds. 

Fig. 8. A snapshot of the x-component of the electric induction, D x for the dielectric 
grating with dielectric constant e = 4 and thickness h = 0.6/im. The results are presented as 
a function of x and z coordinates measured in units of D. Red and blue colors correspond, 
respectively, to positive and negative values of the induction. The snapshot has been produced 
after a sufficiently long propagation time so that the induction pictured in it indeed corresponds 
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to the resonance giving enhanced transmission at A = 1.3D (see Fig. 6). 
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